home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / tmte.cpp < prev    next >
C/C++ Source or Header  |  1995-01-11  |  7KB  |  191 lines
  1.  
  2. //#define WANT_STREAM
  3. #define WANT_MATH
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmatap.h"
  8.  
  9. void Print(const Matrix& X);
  10. void Print(const UpperTriangularMatrix& X);
  11. void Print(const DiagonalMatrix& X);
  12. void Print(const SymmetricMatrix& X);
  13. void Print(const LowerTriangularMatrix& X);
  14.  
  15. void Clean(Matrix&, Real);
  16. void Clean(DiagonalMatrix&, Real);
  17.  
  18.  
  19.  
  20.  
  21. void trymate()
  22. {
  23. //   cout << "\nFourteenth test of Matrix package\n";
  24.    Tracer et("Fourteenth test of Matrix package");
  25.    Exception::PrintTrace(TRUE);
  26.  
  27.    {
  28.       Tracer et1("Stage 1");
  29.       Matrix A(8,5);
  30. #ifndef ATandT
  31.       Real   a[] =   { 22, 10,  2,  3,  7,
  32.                14,  7, 10,  0,  8,
  33.                -1, 13, -1,-11,  3,
  34.                -3, -2, 13, -2,  4,
  35.             9,  8,  1, -2,  4,
  36.             9,  1, -7,  5, -1,
  37.             2, -6,  6,  5,  1,
  38.             4,  5,  0, -2,  2 };
  39. #else
  40.       Real a[40];
  41.       a[ 0]=22; a[ 1]=10; a[ 2]= 2; a[ 3]= 3; a[ 4]= 7;
  42.       a[ 5]=14; a[ 6]= 7; a[ 7]=10; a[ 8]= 0; a[ 9]= 8;
  43.       a[10]=-1; a[11]=13; a[12]=-1; a[13]=-11;a[14]= 3;
  44.       a[15]=-3; a[16]=-2; a[17]=13; a[18]=-2; a[19]= 4;
  45.       a[20]= 9; a[21]= 8; a[22]= 1; a[23]=-2; a[24]= 4;
  46.       a[25]= 9; a[26]= 1; a[27]=-7; a[28]= 5; a[29]=-1;
  47.       a[30]= 2; a[31]=-6; a[32]= 6; a[33]= 5; a[34]= 1;
  48.       a[35]= 4; a[36]= 5; a[37]= 0; a[38]=-2; a[39]= 2;
  49. #endif
  50.       A << a;
  51.       DiagonalMatrix D; Matrix U; Matrix V;
  52. #ifdef ATandT
  53.       int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
  54. #else
  55.       DiagonalMatrix I(A.Ncols());
  56. #endif
  57.       I=1.0;
  58.       SymmetricMatrix S1; S1 << A.t() * A;
  59.       SymmetricMatrix S2; S2 << A * A.t();
  60.       Real zero = 0.0; SVD(A+zero,D,U,V);
  61.       DiagonalMatrix D1; SVD(A,D1); Print(DiagonalMatrix(D-D1));
  62.       Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
  63.       Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
  64.       Matrix B = U * D * V.t() - A; Clean(B,0.000000001);Print(B);
  65.       D1=0.0;  SVD(A,D1,A); Print(Matrix(A-U));
  66.       SortDescending(D);
  67.       D(1) -= sqrt(1248); D(2) -= 20; D(3) -= sqrt(384);
  68.       Clean(D,0.000000001); Print(D);
  69.       Jacobi(S1, D, V);
  70.       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  71.       SortDescending(D); D(1)-=1248; D(2)-=400; D(3)-=384;
  72.       Clean(D,0.000000001); Print(D);
  73.       Jacobi(S2, D, V);
  74.       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  75.       SortDescending(D); D(1)-=1248; D(2)-=400; D(3)-=384;
  76.       Clean(D,0.000000001); Print(D);
  77.  
  78.       EigenValues(S1, D, V);
  79.       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  80.       D(5)-=1248; D(4)-=400; D(3)-=384;
  81.       Clean(D,0.000000001); Print(D);
  82.       EigenValues(S2, D, V);
  83.       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  84.       D(8)-=1248; D(7)-=400; D(6)-=384;
  85.       Clean(D,0.000000001); Print(D);
  86.  
  87.       EigenValues(S1, D);
  88.       D(5)-=1248; D(4)-=400; D(3)-=384;
  89.       Clean(D,0.000000001); Print(D);
  90.       EigenValues(S2, D);
  91.       D(8)-=1248; D(7)-=400; D(6)-=384;
  92.       Clean(D,0.000000001); Print(D);
  93.    }
  94.    {
  95.       Tracer et1("Stage 2");
  96.       Matrix A(20,21);
  97.       for (int i=1; i<=20; i++) for (int j=1; j<=21; j++)
  98.       { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 21-i; else A(i,j) = -1; }
  99.       A = A.t();
  100.       SymmetricMatrix S1; S1 << A.t() * A;
  101.       SymmetricMatrix S2; S2 << A * A.t();
  102.       DiagonalMatrix D; Matrix U; Matrix V;
  103. #ifdef ATandT
  104.       int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
  105. #else
  106.       DiagonalMatrix I(A.Ncols());
  107. #endif
  108.       I=1.0;
  109.       SVD(A,D,U,V);
  110.       Matrix SU = U.t() * U - I;    Clean(SU,0.000000001); Print(SU);
  111.       Matrix SV = V.t() * V - I;    Clean(SV,0.000000001); Print(SV);
  112.       Matrix B = U * D * V.t() - A; Clean(B,0.000000001);  Print(B);
  113.       for (i=1; i<=20; i++)  D(i) -= sqrt((22-i)*(21-i));
  114.       Clean(D,0.000000001); Print(D);
  115.       Jacobi(S1, D, V);
  116.       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  117.       SortDescending(D);
  118.       for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
  119.       Clean(D,0.000000001); Print(D);
  120.       Jacobi(S2, D, V);
  121.       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  122.       SortDescending(D);
  123.       for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
  124.       Clean(D,0.000000001); Print(D);
  125.  
  126.       EigenValues(S1, D, V);
  127.       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  128.       for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
  129.       Clean(D,0.000000001); Print(D);
  130.       EigenValues(S2, D, V);
  131.       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  132.       for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
  133.       Clean(D,0.000000001); Print(D);
  134.  
  135.       EigenValues(S1, D);
  136.       for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
  137.       Clean(D,0.000000001); Print(D);
  138.       EigenValues(S2, D);
  139.       for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
  140.       Clean(D,0.000000001); Print(D);
  141.    }
  142.  
  143.    {
  144.       Tracer et1("Stage 3");
  145.       Matrix A(30,30);
  146.       for (int i=1; i<=30; i++) for (int j=1; j<=30; j++)
  147.       { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 1; else A(i,j) = -1; }
  148.       Real d1 = A.LogDeterminant().Value();
  149.       DiagonalMatrix D; Matrix U; Matrix V;
  150. #ifdef ATandT
  151.       int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
  152. #else
  153.       DiagonalMatrix I(A.Ncols());
  154. #endif
  155.       I=1.0;
  156.       SVD(A,D,U,V);
  157.       Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
  158.       Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
  159.       Real d2 = D.LogDeterminant().Value();
  160.       Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
  161.       SortDescending(D);  // Print(D);
  162.       Real d3 = D.LogDeterminant().Value();
  163.       ColumnVector Test(3);
  164.       Test(1) = d1 - 1; Test(2) = d2 - 1; Test(3) = d3 - 1;
  165.       Clean(Test,0.00000001); Print(Test); // only 8 decimal figures
  166.       A.ReDimension(2,2);
  167.       Real a = 1.5; Real b = 2; Real c = 2 * (a*a + b*b);
  168.       A << a << b << a << b;
  169.       I.ReDimension(2); I=1;
  170.       SVD(A,D,U,V);
  171.       SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
  172.       SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
  173.       B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
  174.       D = D*D; SortDescending(D);
  175.       DiagonalMatrix D50(2); D50 << c << 0; D = D - D50;
  176.       Clean(D,0.000000001);
  177.       Print(D);
  178.       A << a << a << b << b;
  179.       SVD(A,D,U,V);
  180.       SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
  181.       SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
  182.       B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
  183.       D = D*D; SortDescending(D);
  184.       D = D - D50;
  185.       Clean(D,0.000000001);
  186.       Print(D);
  187.    }
  188.  
  189. //   cout << "\nEnd of Fourteenth test\n";
  190. }
  191.